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^1 ' Abstract. We have developed a numerical code with which we study 

Q_|. the effects of 2D perturbations on stellar structure. We present new 

^ | numerical and analytical results on the heating of a main-sequence star 

■* ■ in a binary system by its companion. 



J> ■ 1. Introduction. 

in 

Numerical study of stellar structure and evolution is now a well established 
subject with a long history. Today there is a plethora of ID stellar evolution 
codes, from which astrophysicists can choose. All these models rely on the 
£Sj \ assumption that the star is spherically symmetric. 

■ This assumption has served well for many years. However it is insufficient in 

^ (— I cases that require stars to be non-symmetric. Stellar rotation, accretion and bi- 

Q H i nary systems all add asymmetry to a spherical star. The required computational 

power is now available to look at these perturbations. 
5— i ' We have been working to develop a code that calculates 2D stellar structures 

in a non-computationally demanding way. So far we have studied perturbations 
to structure only but we plan to move into 2D stellar evolution. The two areas we 
have studied are rotation of solar-like stars and heating of a star by a companion 
star. The results are presented in this paper. 
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2. Numerical Method. 

The method we use was developed by Garaud (2002). A ID model in hydrostatic 
equilibrium is taken from a standard ID evolution code and used as the under- 
lying state for perturbation to the structure in 2D. The perturbative equations 
are linear with respect to thermodynamic quantities, with full non-linearity kept 
for the flow perturbations. The equations are 

pu ■ Vu = - Vp - pV<J> - p~V$ + zA7 2 u + -iA7(V • u), (1) 
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V 2 <l = AnGp, (2) 

+ (3) 
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V • (pu) = 0, 



(4) 



and 



pTu • Vs = kV 2 f, 



(5) 



where u is fluid velocity, p the pressure, p the density, $ the gravitational poten- 
tial, v the kinematic viscosity, T the temperature, k the thermal conductivity 
and s the entropy. The variables with a tilde are the perturbed variables and 
those without are taken from the ID hydrostatic equilibrium model. We look 
for solutions to these equations that are stationary and axisymmetric. Note that 
we do not perturb the energy generation or opacities. The solution method is 
to use finite differences in the radial direction and expansion upon Chebyshev 
polynomials (essentially a Fourier decomposition) in the latitudinal direction. 

Once this system is setup all that remains is to impose the boundary con- 
ditions appropriate to the system considered. 

3. Stellar Rotation. 

Stellar Rotation is an old problem (Eddington 1925; Sweet 1950). Commonly 
used models examine meridional currents in stars that are rotating slowly, that 
is more slowly than the Sun. There are two main effects of rotation on stellar 
evolution, the increase or decrease of the amount of contraction or expansion of 
a star during its evolution and the induction of extra mixing of material within 
the star, especially within regions that are convectively stable. This has been 
recently reviewed by Pinsonneault (1997). 

Garaud (2002) studied the meridional currents in the radiative zone of solar- 
like stars set up by the differential rotation of an outer convective zone. To 
investigate the problem a radiative zone model was taken from a ID standard 
solar model (Christensen-Dalsgaard et. al. 1991) and perturbed by forcing ro- 
tation on the outer boundary and the equations were solved within the radiative 
zone of the star, with the assumption that the convection is not affected by the 
radiative zone because the driving forces are larger than any effect from the 
radiative-zone flow. 

The boundary condition imposed differential rotation at the surface of the 
radiative zone of the form 



along with stress continuity across the surface and V 2 T = outside. For further 
details see Garaud (2002). 

4. Heating of Main- Sequence Stars in X-ray Binaries. 

It is possible to study different problems with the same method. Essentially 
only the boundary conditions require alteration. Therefore we adapted the code 
to study heating of a main-sequence star in a binary system by X-rays emitted 
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by an accreting neutron-star companion. Many such X-ray binaries are known 
to exist. 

Observational evidence that something non-ID is occurring comes from 
HZ Herculis where there are variations in the light curve, especially at the 
minimum when the X-ray source is eclipsed, that cannot be explained by a 
main-sequence star with a homogeneous surface temperature. It is generally 
assumed that the external radiation that falls on the star's surface is thermal- 
ized and re-emitted, with some of the energy emitted from the unilluminated 
side. To achieve this there must be some form of circulation current in the star's 
envelope. 

There have been previous attempts to model this situation, Kippenhahn & 
Thomas (1978) assumed that the flux was thermalized very near to the surface 
and found that the induced flow penetrated to a depth of 20% of the star's 
radius, with a 10% increase in the illuminated hemisphere temperature and 2% 
of the extra flux emitted from the unilluminated surface. The velocities of the 
vertical flows are subsonic, while the flows parallel to the surface are likely to 
be supersonic. 

Later, Kirbiyik (1982) produced a similar model but with a much shallower 
penetration of only a few percent of the star's radius. These analyses of the 
reflection effect have been used to explain many light curves in close binary 
systems, with the solution of Kippenhahn & Thomas (1978) being the most 
commonly citied. 

4.1. Boundary Conditions 

Careful thought must be given to formulation of the boundary conditions because 
there are many possibilities for boundary conditions that add extra energy to 
the surface of the star. The X-ray flux is absorbed in a thin layer near the 
surface. The column density required is of the order of 0.1 — > lgcm -2 . This 
is reached at a very shallow depth in the stellar surface, under the assumption 
that there is no absorption in the stellar atmosphere. 

To a first approximation we can say that, at the surface, we are adding 
extra energy with a boundary condition 

cfT 

-k— = a B T 4 -af (9) cosO. (7) 
The function / is defined as 

f t Q\ f 1 ~ 2 ^ ^ ^ 2> (S\ 

J ^ \ otherwise, ^ ' 

where 9 is the latitude, a the strength of the X-ray flux, a = Lx/i^d 2 ). 
Because we are dealing with small perturbations we expand T — > T + T and 
assuming T is small we obtain 

df 

-k— = Aa B T 3 f - af(6) cos 9. (9) 

This boundary condition can also be used if the illuminating star is a main- 
sequence star. 
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We also have conditions on the flow such that 

u r (r = 1) = and u e {jx = ±1) = 0, (10) 

along with a stress free boundary on the flow at the surface. We have defined 
H = cos 0, r the dimensionless radial coordinate, r = is at the centre and r = 1 
is at the surface. 

4.2. Numerical Results 



"i 1 r 

Numeric 1st Order - 
Numeric 2nd Order - 
Analytic 1st Order - 
Analytic 2nd Order - 



-10 

0.86 




0.88 0.9 0.92 0.94 
Fractional Radius 



0.96 



0.98 



Figure 1. Plot of the radial dependence of the temperature pertur- 
bations for the solar radiative zone (surface taken to be at 0.6Rq), 
the illuminating X-ray flux is 10 35 ergs _1 at 1 au. Notice the shallow 
penetration and the good agreement of analytic and numerical results. 



Only one model, the solar radiative zone, has been studied successfully 
numerically. By this we mean that we have removed the upper part of the 
Sun so we are directly heating the region at 0.6Rq. This was done to test the 
numerics on a well-known background state and helped us develop analytical 
tools that can be generalised to any stellar radiative zone. The viscosity has 
also been boosted to resolve the thin thermo-viscous boundary layers. 

We find that circulation is induced but limited to a thin boundary layer 
near the surface. Similarly to previous authors we find that there are two cir- 
culation cells. Figure 1 shows the temperature fluctuations versus depth. The 
perturbations only penetrate to a shallow depth. Figure 2 plots the value of ip 
which represents the nature of the induced flows. Radial velocities are subsonic 
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but latitudinal velocities do reach supersonic values for illumination equivalent 
to a source of 10 35 ergs _1 at 1 au. The flow velocities are linearly dependent on 
the illuminating flux. 

Because the effect of the heating is confined to a small boundary layer we 
realised that we could look at the problem analytically. This provides us with 
some generic analytical results on the heating that we can use for any star and 
any form of heating. 

4.3. Analytic Modelling 

We have investigated the situation and solved the equations analytically using 
boundary layer approximations to keep only the highest order radial derivatives. 
We use the same method of Fourier expansion of the variables. We simplify the 
equations by restricting them to the surface layer and using only the dominant 
terms in that region. We start by solving the mass conservation equation for 
incompressible flow 

V • (pu) = 0. (11) 

We assume the density variation is negligible and define a circulation variable tp 
such that 

1 dip ^ 1 dip 

H dfi ry/1 - [i 2 or 

We now make the variables dimensionless by scaling distance to the star's radius, 
R, and velocities using the viscous diffusion timescale so that tp' = ipRu becomes 
dimensionless. Thus the energy equation (5) and vorticity equation (1) become 

(l-U 2 )^--^- (13) 

and 

uRpc p N 2 Tdip' o~ , 

— ^ = V 2 T. (14) 

kg o\x 

where g is the acceleration due to gravity, c p the specific heat capacity at constant 
pressure and iV 2 the buoyancy frequency. 

To solve the equations, we expand the variables T and ip as 

f - T {r) + /iTi(r) + (2fi 2 - l)T 2 (r) (15) 

and 

«o(r) + + (2/T - l)u 2 (r). (16) 



To zeroth order the system reduces to Laplace's equation for the tempera- 
ture perturbation. For the higher orders we use equations (12)-(16) to obtain 
equations of the form 

Ti = ^ (17) 
when % > 0. These have a solution of the form 



Ttexpfar-lj), (18) 
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Latitude. 



Figure 2. Plots of T(r, 6) and ip for the solar radiative zone model 
and the numerical method, the details are identical to those in figure 1. 
Greyscale levels represent T(r,6). With white the highest temperature 
and black the lowest. Line contours of ip represent the circulations. 
Solid lines are clockwise flows, dashed anti-clockwise. The velocities 
are greatest at the surface and decay with depth as determined by 5. 
The box is 2.6 x 10 11 by 6.5 x 10 9 cm. Scaled to real viscosity the y-axis 
goes from 1.0 to 0.99990. 



with similar equations for the flow velocities. The value of 8 is of most interest 
because it determines the depth of penetration below the surface. With this 
analysis we obtain 

'.-( 3a ^)'-*-( SH ^)*- 

This gives us six solutions for Si in the complex plane and we use those that 
decay with depth. We find the full solution becomes 

Ti(r) = T^exp far - 1)) + 2exp (^(r - 1)) cos (^(r - 1))^. (20) 

Then by fitting to the boundary condition (9) we obtain 

^ = A l a(l2a B T 3 + ^)" 1 . (21) 

Orders only differ by a numerical factor, for expansion to second order 
A\ = \ and A2 = The solutions for the flow components u r and uq can be 
deduced from equations (13) and (21). Table 1 shows values of S\ for various 
stellar models, computed with the Eggleton code (see Pols et al. 1996). The stars 
are above 2Mq to ensure fully radiative envelopes. In the numerical study we 
had to boost the viscosity by a few orders of magnitude to resolve the boundary 
layer. But in the analytic work we use real values for the viscosity (Spitzer 1962) 
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and so the penetration is shallower, varying as 5 oc However if we use the 
same values as for the numerical calculations the results agree closely. See figure 
1 for the comparison of results by this method and those from numerical results. 

The temperature perturbations are also listed in table 1, with example flow 
velocities in table 2. The sources are placed at a distance of 1 au from the 
secondary. This is much larger than the typical separation in an X-ray binary (a 
few Rq but we cannot yet deal with such intense irradiation - see Discussion) . 
For one of the stellar models we show the variation of the values with different 
illuminating fluxes. 



Table 1. Typical values of 5\ and % from the analytic method. Stel- 
lar models were obtained from the Eggleton stellar evolution code. 5\ 
is dimensionless, R5± 1 is the penetration depth. 



Star Mass 




R 


Penetration 


Lx 


% 


Ti 


T 2 






/% 


depth / Rq 


/ergs -1 


/mK 


/mK 


/mK 


2M 


2962 


2.339 


0.00079 


10 33 


0.671 


0.193 


0.0750 








10 35 


67.1 


19.3 


7.50 










10 37 


6710 


1930 


750 


3M 


1965 


3.001 


0.00153 


10 35 


25.1 


6.59 


2.54 


4M 


1866 


3.417 


0.00183 


10 35 


15.8 


4.12 


1.59 


6M 


1790 


4.434 


0.00248 


10 35 


6.73 


1.66 


0.63 


8M 


1719 


5.055 


0.00294 


10 35 


4.55 


1.11 


0.43 


1OM 


1657 


5.769 


0.00348 


10 35 


3.16 


0.76 


0.29 



Table 2. Maximum values of velocities, in terms of the sound speed, 
c s . Values are from the surface, velocity decays with depth. 



Star Mass 


Cg/ kms 1 


Lx/ergs~ x 


Ug(fJ, = 0.3)/c s 


2M 


12.8 


10 33 


2.99 x 10~ 2 




10 35 


2.99 






10 37 


299 


3M 


16.5 


10 35 


1.211 


4M 


17.6 


10 35 


0.780 


6M 


21.4 


10 35 


0.225 


8M 


23.0 


10 35 


0.136 


1OM 


24.3 


10 35 


0.085 



4.4. Discussion 

From the tables it can be seen that a circulation is present even if the irradiating 
flux is low, although the perturbation is too small to be observed. With larger 
fluxes a problem arises because the latitudinal flow becomes supersonic. This 
agrees with the results of Kippenhahn and Thomas but would produce shocks 
in the flow for which we have not allowed in this work. 
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The amount of extra flux emitted from the unilluminated side of the star 
can be estimated by comparison of the strength of the dipole mode to the other 
terms. We find that about 60% of the flux is emitted from the irradiated surface 
with 40% from the unilluminated side. This is a smaller difference than found by 
previous authors and is due to the large zeroth order temperature perturbation. 
The value is large because we have not allowed the star to expand. 

There are two interesting deviations from the work of Kippenhahn and 
Thomas. First we have much shallower penetration similar to the results of 
Kirbiyik and second there is more flux emitted from the opposite side. In our 
next step we shall test these analytic results for various stars numerically. 

In these results it should also be noted that the boundary conditions are 
simple. While they are a reasonable estimate for illumination by black-body-like 
stars, if the bulk of the incident radiation is in the form of X-rays then energy 
will not all be deposited at the surface of the star but at a greater depth. 

5. Future Work 

The code calculates 2D perturbations to the hydrostatic equilibrium of any stel- 
lar radiative region. However it can only deal with linear perturbations to the 
thermodynamic variables. Therefore the next phase is to develop the code to 
model non-linear effects and solve for the full structure, not just perturbations 
to an assumed hydrostatic equilibrium. In particular in the above study, we have 
not allowed the star to expand or shrink when heated and we shall need to allow 
deformation, deviations from spherical shape. We shall also improve upon the 
boundary conditions to include a model atmosphere as in Tout et al. (1989) and 
to include external radiation pressure from the irradiating flux, along with the 
possible deeper penetration of the incident radiation. To study less massive stars 
we shall also be looking at expanding the code to deal with convection zones by 
some parameterisation that does not require us to follow the full dynamics of 
convection. 
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